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1 ABSTRACT 

ph ' Using time evolutions of the relevant linearised equations we study non-axisymmetric 

O |. oscillations of rapidly rotating and superfluid neutron stars. We consider perturbations 

of Newtonian axisymmetric background configurations and account for the presence 
of superfluid components via the standard two-fluid model. Within the Cowling ap- 
^ ' proximation, we are able to carry out evolutions for uniformly rotating stars up to 

. the mass-shedding limit. This leads to the first detailed analysis of superfluid neutron 

star oscillations in the fast rotation regime, where the star is significantly deformed 
, by the centrifugal force. For simplicity, we focus on background models where the two 

J> ' fluids (superfluid neutrons and protons) co-rotate, are in /3-equilibrium and coexist 

0^ , throughout the volume of the star. We construct sequences of rotating stars for two 

analytical model equations of state. These models represent relatively simple general- 
isations of single fluid, polytropic stars. We study the effects of entrainment, rotation 
and symmetry energy on non-radial oscillations of these models. Our results show 
that entrainment and symmetry energy can have a significant effect on the rotational 
qq ' splitting of non-axisymmetric modes. In particular, the symmetry energy modifies the 

, inertial mode frequencies considerably in the regime of fast rotation. 



o 



in 

m 



x 



Key words: methods: numerical - stars: neutron stars: oscillation star:rotation. 



1 INTRODUCTION 

According to the standard paradigm, millisecond pulsars are accelerated to their fast rotation rates by accreting matter 
from a close companion. This means that they tend to be relatively old. Moreover, the fastest spinning pulsars should have 
weak (exterior) magnetic fields. In the standard accretion model, neutron stars with canonical 10 12 G dipole fields will reach 
equilibrium already at a modest spin. A weak surface field is also expected since accretion leads to magnetic field burial. This 
picture agrees well with observational data. Rapidly rotating neutron stars are most commonly found in binary systems. It 
is well established that accreting neutron stars in low-mass X-ray binaries (where the angular momentum transfer is more 
efficient due to the long evolution time of the low mass partner) can r each a millisecond r otation period. Furthermore, the 
fastest known millisecond pulsar J1748-2446ad, with a period of 1.39 ms (He ssels et alj |2006h is in a binary system. In fact, its 



companion, with mass M O.14M0, could still fill the Roche lobe powering the spin-up phase further. The mass and radius 
of J1748-2446ad are unknown, but combining reasonable ranges for these parameters, M = 1.4 — 2A4q and R = 10 — 14 km, 



with an empirical formula for the maximum rotation of the star (jLattimer fc Prak ash 2004) 



M \ 1/2 / 10km \ 3/2 



Qk « 6566 ^— J (-g-j Hz, (1) 

one finds that the spin of this object lies in the range 0.48 < Q/Qk < 0.96. In other words, it could be close to the 
mass-shedding limit. 
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The temperature of a mature neutron star is likely below the critical temperature, T c ~ 5 x 10 9 K, where neutrons 
and protons are superfluid and superconducting, respectively. Depending on the cooling mechanism, neutrino emission can 
cool a hot proto-neutron sta r below this temperature shortly after its formation in a core collapse supernova (see e.g, 
Page. Geppert fc WebeJbood ). In an accreting system, the neutron star core temperature is not expected to increase be- 



yond T c (nuclear burning in the accreted surface layers is thought to the heat the core to ~ 10 8 K) . Hence, all mature neutron 
stars should contain degenerate superfluid neutrons in the outer core and the inner crust and degenerate superconducting 
protons in the outer core. The deep core may contain more exotic phases of matter, like superfluid hyperons and/or colour 
superconducting quarks. Superfluidity influences the thermal evolution and the dynamical properties of a neutron star. In 
particular, the dynamics is strongly affected by entrainment, the formation of quantized neutron vortices, and the presence 
of new dissipative mechanisms like mutual friction. An understanding of superfluid dynamics is crucial for modelling many 
aspects of the neutron star physics, ranging from pulsar glitches and free precession to the mutual friction damping of stellar 
oscillations and associated instabilities. 



Tidal forces, accretion and glitches may trigger oscillations and/or instabilities in rapidly rotating neutron stars. Obser- 
vations of such oscillati ons, either via electromagnet i c or gravitational radiation, would help us explore the exotic physics of 
these compact objects ( Andersson fc Kokkotasl 19981 : Benhar et al. 20041 : Samuelsson fc Andersson 2007 ). In this context it is 
interesting to note the differences in dynamics between neutron stars above the superfluid transition temperature and colder 
systems. The core of a hot young neutron star is relatively well described by the Navier-Stokes equations. In contrast, a star 
with a superfluid core requir es a multi-fluid description. The standard model for such system s is inspired by the two-fluid 
model for superfluid Helium ijLandau fc Lifshita Il959l : lKhalatnikovlll965l : iTillev fc Tillevlll990h . The oscillation spectrum of 



cv 

superfluid neutron stars reflects the presence of the additional degree of freedom ( Epstein! 198a ). Basically, the perturbed fluid 
elements i n a two- fluid system can osci ll ate either in phas e or in a counter-phase motion . Previous studies, see for example, 



Leel (|l995T ); Undersson fc Comerl |200ll ): IPrix fc Rieutordl |2002h : lYoshida fc Led (|2003al ) . have established that co-moving 
pulsations have spectral properties similar to single fluid stars. Hence, it is natural to refer to such modes as "ordinary" 
modes. The counter-moving degree of freedom leads to new oscillation modes that are specific to the two-fluid systems. These 
are often referred to as "superfluid" modes. Later, when we write down the perturbation equations for a superfluid neutron 
star core, we choose to work with variables that are directly associated with the two types of motion. This is natural if we 
want to distinguish spectral properties associated with the "superfluid" degree of freedom. In addition, we use the standard 
classification of neutron star oscillation modes, based on the main restoring force that acts on a perturbed fluid element. 
A rotating single fluid star can sustain acoustic and inertial modes restored by pressure and the Coriolis force, respectively. 
When t hermal or compo s ition gradients are present in the star, buoyancy acts as restoring force for the so-called gravity 
modes ( Unno et al ] |l989l : iReiseneeeer fc Goldreichlll992l ). Previous work has shown that superfluid neutron stars have two 
families of acoustic and inertial modes, more or less clearly (depending on the stellar model) associated with the co- and 
counter-moving fluid motion. It is, however, not the case that al l single flu i d modes have a "double" in the superfluid problem . 
The gravity modes are not present at all in a superfluid core (|Ledll995l ; I Andersson fc Comerl 1200 ll ; IPrix fc Rieutordl 120021 ). 
Their absence provides a potentially important signature for neutron star seismology. 



Rapidly rotating neutron stars have been studied in detail with a variety of methods (see e.g., IStergioulasll2003l ). Yet, 
there have not been any previous studies of multi-fluid dynamics in the rapid rotation regime near the break-up limit. The 
oscillation mod es of superfluid neutron st a rs have only been cal culated in the frequency domain using the slow rotation 
approximation ( Lindblom fc Mendelll 2000l : Yoshida fc Lee! 2003al lbh . In that framework, the effects of the stellar rotation 
are determined perturbatively as corrections to the non-rotating results. T he only previous attempt (t hat we are aware of) 
to study superfluid oscillat ions for truly fast spinni ng stars is the work bv lLindblom fc Mendelll (2000). They extended the 
two-potential formalism of Ipser fc Lindblom 1 1990h to the superfluid case. However, due to numerical difficulties they could 
not study rotating models near the break-up limit. Neither did they manage to determine the superfluid modes. 



In this work, we study the time evolution of perturbed fast rotating, Newtonian superfluid neutron stars within the 
Cowling approximation. As far as we are aware, this is the first study that evolves in time the oscillations of superfluid 
neutron stars. Moreover, it is the first detailed analysis of the rapid rotation regime. Within the framework of the two- 
fluid formalism, we carry out a linear perturbation analysis for stationary and axisymmetric equilibrium configurations. As 
preparation for more detailed studies, we consider relatively simple models where the two fluids coexist throughout the star, 
and where the unperturbed configuration is in /3-equilibrium. These assumptions imply that the two fluids are co-rotating 
and share the same external stellar surface in our background configurations. We use these models to investigate the effects 
of entrainment, symmetry energy and rotation on the superfluid oscillation spectrum. In order to establish the reliability of 
our numerical evolution code, we compare our results to previous work for non-rotating and slowly rotating models. We then 
consider, for the first time, the dynamics of superfluid models rotating up to the mass shedding limit. 
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2 EQUATIONS OF MOTION 

We use the two-fluid framework for superfluid neutron stars (Mendel] 1991al lbl; IPrix 2004 ; Andersson fc Comer 200(j) . This 
model distinguishes between a superfluid neutron component and a neutral mixture of protons and electrons. The charged 
particles are assumed to be locked together by the electromagnetic interaction on a timescale much shorter than the dynamics 
we consider. For simplicity, we refer to the charged particle conglomerate as the "protons" from now on. When the mass of 
each component is conserved the dynamics is d escribed b y two mass conservation laws, two Euler-type equations and the 
Poisson equation for the gravitational potential (jPrixl 120041 ): 



<9t/O x + V; (/5 X 



0. 



ft 
Px 



(2) 



(3) 



V 2 $ = 47rGp. (4) 

Here the indices i and k label the spatial components of the vectors while x and y denote the two fluid components. The 
constituent indices take the values n for the neutrons and p for the protons. Note that the summation rule for repeated indices 
applies only for spatial and not for constituent indices. In equations ©-Q, p = p n + pp and $ represent the total mass 
density and the gravitational potential, respectively. Meanwhile, f x is the force density acting on the x fluid component. In 
this work, we consider an isolated system where dissipation processes, like mutual friction, are absent. We then have f x = 0. 
Furthermore, we have assumed that the particle masses are equal, m = m n = m p , and defined the chemical potential and the 
relative velocity between the two fluids as follows: 

d£ 

''/'•: „. W 2 ' 



/<x 



(5) 
(6) 



The energy functional £ = £(n n ,n p ,w np ) describes the equation of state (EoS) of the system. Finally, the non-dissipative 
entrainment between the two fluids is governed by the parameter e x , which follows from the definition: 

d£ 



2p x 



dwl p 



(7) 
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The equations that describe rapidly and uni formly rotating background models can be derived by integrating the Euler- 
type equations ([3]) and the Poisson equation (|4]) ijPrix et al.ll2002l : lYoshida fc Eriguchill2o"oi ). This leads to 



/ix 



(8) 



where fl x and C x are, respectively, the angular velocities and the integration constants for the neutron and proton fluids. In 
this work, we focus on corotating background models, Q n — Q p — fl, where the two fluids are in /3-equilibrium and have a 
common surface. The hydrostatic equilibrium equation |S} then become 

£ + — sin(9 2 ft 2 = C, (9) 



where p, p = fia = p, is the background chemical potential and C = C n — C p . It is worth noticing that equations Q 
and j9l) are formally equivalent to the single fluid problem provided one replaces the chemical potential with the en- 
thalpy jYoshida fc Erigu chi 2004) . Giv e n an equation of state, equations Q and ((9} can be numerically solved using the 
self-consistent f ield method of^ HachisuJ ( 19861 ) . The surface of the star corresponds to the zero chemical potential surface, 
ft (r (0) , 0) = (jYoshida fc Eriguchilbool T" 



3 PERTURBATION EQUATIONS 

It is straighforward to write down the system of partial differential equations that governs the Eulerian perturbations 
#p x , fa x, Stix and However, instead of working with these variables, we define (see lAndersson. Glampedakis fc Haskell 



2008, for a detailed discussion) new variables which are more closely related to the natural degrees of freedom of the problem. 



In the rotating frame, the co-moving (ordinary) motion is described by the mass flux f (not to be confused with the mutual 
force f x ), the total mass density Sp and the pressure SP, defined by 

f = ppVp + p n v n , (10) 
Sp = Sp n + Spp, (11) 

VSP = Pp VSp p + p n VSp n . (12) 
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Meanwhile, the counter-moving (superfluid) motion is described by a vector field D that is proportional to the relative velocity 
between the two fluids, the scalar perturbation 6/3 that measures the deviation from /^-equilibrium and the quantity 8x P , which 
is related to the perturbed proton fraction. These variables are defined by 

D = x p (1 - x p ) p (v p - v n ) , (13) 
5(3 = 8p p - 6p n , (14) 
SXp = pSx p , (15) 

where x p — p p /p is the proton fraction. 

In order to simplify the evolutions, we neglect the perturbations of the gravitational potential <5<!> = 0. That is, we adopt 
the Cowling approximation. This approximation should be quite accurate for inertial modes. For low order acoustic modes, 
like the f-mode, it is not so accurate but the results are still qualitatively correct. For our present purposes, this should be 
sufficent. Although, it would not be too difficult to solve also for the perturbed gravitational potential it is computationally 
costly to add the solution of an elliptic equation to our evolutions. Hence, we decided not to solve the full problem in this 
first exploratory study. 

In the frame of the rotating background the final perturbation equations can be written 

VP 

d t f = —\7SP — 2fi x f + — — Sp , (16) 

P 

x p {\-x p )pVS(3 2tixD 

9tD = i=i r~r ' (17) 

d t S P = -V-f, (18) 
d t S Xp = -V-D-f-V:r p , (19) 

where we have defined e = e n /x p . 

The time evolution of the non-axisymmetric perturbation equations is a three-dimensional problem in space. However, 
linear perturbations on an axisymmetric b ackground can be expanded in terms of a set of basis functions (cos m<j) , sin mcj)), 



where m is the azimuthal harmonic index i Papaloizou fc Pringkj|l98C ). The mass density perturbations as well as the other 
perturbation quantities take the following form rtjones et alj|2002l ; IPassamonti et al. 20091 ) 



5p(t,r,6,<fi) = \5ptn (t, r i 9) cosm<j> + Sp m (t, r, 6) sinm^] . (20) 

m = 

The perturbation equations now decouple with respect to m and the problem becomes two-dimensional. Therefore for any m, 
the non-axisymmetric oscillations of a superfluid neutron star requires the solution of a system of eighteen partial differential 
equations for the twenty variables (f , 5p , SP , D , 5xp,8p ). To fully specify the problem, the set of equations (|16|l - (|19[l 
must be complemented by two relations that depend on the EoS (see Section [4|. 



3.1 Boundary Conditions 

In order to evolve the perturbation equations we must also specify boundary conditions. For non-axisymmetric oscillations 
with m ^ 2, equations (|16p - (|19|) are regular at the origin, r — 0, when the following conditions are satisfied: 

5P = 8 Xp = 5/3 = 5p = , and f = D = . (21) 

For the boundary condition at the stellar surface, it is worth remembering that the unperturbed configuration is such that 
the two fluids have a common surface (see Section [TJ. At the perturbed level, we require that the Lagrangian perturbation of 
the chemical potentials vanishes at the surface, i.e., 

A x /i x = Sp* + £ x ■ V/i x = , (22) 

where the Lagrangian variations A x are associated with the fluid displacements £ x ([Andersson. Comer fc Grosartl |200J) . 
Equations (|22[l can be expressed in terms of SP and S/3 by using the definitions (I12|l and (|14[l . This leads to 

AP = *P+[ipCp + (l-Xp)ej- VP = 0, (23) 

A/3 = 60 + (Z v -Zn)- — =O, (24) 

P 

where we have used the fact that the background model is in /3-equilibrium condition, i.e., p n = p P - As in the single fluid 
case we can derive a simpler condition for the pressure perturbation. From the Euler equation for the stationary background 
(noting that (j 12 1) holds also at the unperturbed level) 

VP = -pfl x (O x r) - pV$, (25) 
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and the fact that the EoS us ed in this wor k are such that the total mass density vanishes at the surface, it follows that 
equation ((23]) is equivalent to ( Tassoul 19781 ): 

SP = . (26) 
This condition is numerically convenient, and ensures that all variables remain regular at the surface. 



The reflection symmetry with respect to the equator divides the perturbation variables into two sets (|Passamonti et al 



2009). Perturbations of the Type I parity class have (/ r , Sp, SP, D r , D^, <5x P , Sf3) even and (fe,Dg) odd. The opposite is 



true for perturbations of Type II. 



4 EQUATION OF STATE 

To complete the formulation of the superfluid oscillation problem we need to supply a suitable multi-parameter equation 
of state. In a truly realistic model, the equation of state should be obtained from a microscopic (quantum) analysis. It will 
completely specify not only the relation between pressure and density, but also the composition and the detailed superfluid 
energy gaps for neutrons and protons. Moreover, it has to also provide the entrainme nt parameters. We do not yet have such 
a model, although recent developments in this direction are promising I ChameJ 20081 ). Given this, and the fact that our main 
aim is to understand the oscillations of a rotating superfluid neutron star at the qualitative level, we will opt to work with two 
simple analytic model equations of state. These models, described in Sections 15. H and 15.21 below, are natural generalisations 
of the single fluid polytropes. The analytical models are particularly useful since they allow us to tune key parameters like 
entrainment and symmetry energy more or less freely. As we will see, we can also vary the coupling between the co- and 
counter-moving fluid degrees of freedom. 

Quite generally, the required energy functional can be expressed as a function of the two fluid mass densities and the 
relative velocity: 

£ = £ (Pn,Pp,Wnp) , (27) 

where the dependence on w np ensures Galilean invariance. For a small relative velocity between the two fluids, equation (|27|) 
can be written 

S = So (p n , p P ) + a (p n , P P ) wlp + O (w np ) , (28) 

in which case the bulk EoS So and the entrainment ao can be independently specified at w np = 0. From equation J7]) follows 
that the entrainment parameter £ x is related to the function ao by 

p x£x = 2a . (29) 

Having specified the EoS, we must determine two equations that close the system (|16[) - (|19[) . One possible choice is to 
determine the pressure perturbation and the quantity S(3 from the total density and the proton fraction 

SP = 5P(Sp,Sxp) , 6/3 = 5/3(6p,6x p ) . (30) 

The required relations are obtained by first expressing these quantities in terms of the chemical potential perturbations <5/i x , 
i.e. using the thermodynamic definitions 

SP = PpSflp + pnSpu , (31) 

8/3 = Sjxp-Sfln. (32) 
For corotating background models (w np = 0), the perturbation of the chemical potential p x = /2 X (p p ,p n ) can be expressed as 



dpp 



Spp + Sp n , (33) 



where the mass densities of the two fluid components are defined in terms of the total mass density and proton fraction by 

Spp = XpSp + pSxp , (34) 

8p n = (1 — x p ) Sp — pSxp . (35) 
By introducing equations (|33|) — (|35|) into equations (|31[) - (|32[) . we obtain: 

SP = [((l + 2a)x'p-2(l + a)xp + l) A nn +x 2 p App] pSp+[((l + 2a)xp-l-a)A nn + XpApp]p8xp, (36) 

8(3 = [((1 + 2a) Xp - 1 - <r) A^ + x p A pp ] Sp + [(1 + 2a) A m + A pp ] 5% P , (37) 

where the matrix A xy is defined by 

_ <9/i x _ d 2 S 



dp y dp y dp x 



(38) 
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Table 1. This tabic provides the main parameters for the rotating models A and B (see Sections 15.11 and 15.21 for the detailed EoS). 
The first column labels each model. In the second and third columns we show, respectively, the ratio of polar to equatorial axes and 
the angular velocity of the star. In the fourth column, the rotation rate is compared to the Kepler velocity Qk that represents the mass 
shedding limit. The ratio between the rotational kinetic energy and gravitational potential energy T/|W| and the stellar mass are given 
in the fifth and sixth columns, respectively. All quantities are given in dimensionless units, where G is the gravitational constant, p c 
represents the central mass density and R eq is the equatorial radius. 



Model 


Rp j Req 




n/n K 


T/\W\ x 10 2 


M/(p c R 3 ea ) 


AO 


1.000 


0.000 


0.000 


0.000 


1.273 


Al 


0.996 


0.084 


0.116 


0.096 


1.270 


A2 


0.983 


0.167 


0.230 


0.385 


1.248 


A3 


0.950 


0.287 


0.396 


1.169 


1.197 


Al 


0.900 


0.403 


0.556 


2.385 


1.118 


A5 


0.850 


0.488 


0.673 


3.639 


1.038 


A6 


0.800 


0.556 


0.767 


4.933 


0.956 


A7 


0.750 


0.612 


0.844 


6.252 


0.869 


1 io 


U. 1 uu 


U.UJO 


n Qfi7 




n 77Q 


A9 


0.650 


0.693 


0.956 


8.822 


0.684 


A10 


0.600 


0.717 


0.989 


9.865 


0.579 


All 


0.558 


0.725 


0.999 


10.277 


0.480 


BO 


1.000 


0.000 


0.000 


0.000 


1.833 


Bl 


0.996 


0.094 


0.107 


0.105 


1.825 


B2 


0.983 


0.187 


0.213 


0.419 


1.799 


B3 


0.950 


0.323 


0.368 


1.275 


1.733 


B4 


0.900 


0.453 


0.516 


2.608 


1.632 


B5 


0.850 


0.551 


0.628 


4.002 


1.529 


B6 


0.800 


0.629 


0.717 


5.459 


1.424 


B7 


0.750 


0.695 


0.792 


6.981 


1.317 


B8 


0.700 


0.751 


0.856 


8.566 


1.207 


B9 


0.650 


0.798 


0.909 


11.020 


1.093 


BIO 


0.600 


0.835 


0.952 


11.862 


0.973 


Bll 


0.563 


0.857 


0.977 


13.074 


0.875 


B12 


0.496 


0.877 


0.999 


14.613 


0.669 



and cr corresponds to the so-called "symmetry energy" (|Prix, Comer fc Anderssonl [2003; iPrakash, Lattimer & Ainsworth 

Ann 



19881). That is, we have 



A n 



(39) 



5 STELLAR MODELS 



Background stellar models such that the two fluids are co-rotating can be constructed by solving the hydrostatic equilibrium 
equations Q and for a given bulk EoS £o, cf. (28) . Since equations @ and (|9"| are formall y equivalent to the equilibrium 
equations for a single fluid polytrope, we can straightforwardly use the method of lHachisu (1986) to determine such background 
models (see Section [2). 

For a co-rotating background, entrainment does not not affect the equilibrium configuration. Hence, it can be chosen 
independently from the bulk EoS (see equations 1)280 an d UU)- In fact, the entrainment parameter appears only in the 
perturbation equation (|17l) through the combination 

S = — = e p + e n . (40) 



In the last step we have used the relation (|29p . i.e. p n £n = p P £ P . N uclear physics calculations limit the value of the entrainment 
in the neutron star core to 0.2 ^ £ P ^ 0.8, see Chamell 1 2008! ') for a recent analysis. However, values outside this range 
are pos sible, especial ly for the crust superfluid. In fact, the parameter £ p is expected to have negative values in the crust 
region 1 Chamel 20061 ). Since we are interested in exploring the effect that the different parameters have on the neutron star 
oscillation modes, we will consider the range —0.7 ^ £ ^ 0.7. 
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5.1 Model A 



As our first model equation of state we consider (see, iPrix et alj|2002l ; lYoshida fc Eriguchii 2004) 



(41) 



We will refer to this as model A. Despite its obvious simplicity, this model allows us to investigate the effect of the symmetry 
energy on the oscillation spectrum. This is apparent from the matrix coefficients A xy which have the following form: 

2K 



^4nn 



l-{l + a)x v ' 
2K [1 + a - (1 + 2a) x p ] 
x p [1 — (1 + a)x p ] 
aA nn , 



(42) 

(43) 
(44) 



where a and K are constants. Equations (|42[1 - (|44[) are equivalent to the set used bv lYoshida fc Eriguchii (|2004l ) with 2K = 1/k. 
From the /3-equilibrium condition for the stellar background, jl n = jl p , and the definition (JS| we can derive the proton fraction 
for these models 

(1 + a) A n n 



(l + 2a)A n 



(45) 



If we take the coefficients A xy to be constant, model A leads to a family of non-stratified stars. The perturbation variables 
are related by equations (|36[) - (|37|l which in this case become 

2K(l + a) 



SP = 2KpSp, 



-,$Xp ■ 



(46) 



x p [1 - (1 + a) Xp] 

Therefore, the models are specified by the three parameters K, a and x p . For this class of models, the degrees of freedom that 
describe the co- and counter-moving fluid motion are decoupled. The variables f , 8p, SP evolve independently from D, 8\p, 5/3. 
This is a useful simplification that helps the interpretation of the oscillation spectrum. 

For a background star in /3-equilibrium, the chemical potential is related to the total mass density by 

p = 2Kp, (47) 

which means that the pressure is that of the usual N = 1 polytrope: 

P = Kp 2 . (48) 

We have constructed a family of rotating stars for this equation of state. In dimensionless units, the constant K is 
determi ned automatically by specifying the poly tropic index and the axis ratio between the polar and equatorial radius 



Rp/Req ( Jones et al.ll2002l ; Passamonti et al. 20091) . The properties of our rotating models are given in TablsQl where we also 
label each member of the sequence. The non-rotating model is refered to as AO and the fastest rotating model is All. Being 
non-stratified stars, the values of both the proton fraction x p and the symmetry energy term a affect only the counter-moving 
motion of the perturbed fluid. In fact, they appear only in the perturbation equations (|17p and (|19[) . and in equation (|46l) 
for the 5/3 perturbation. In our numerical simulations, we focus mainly on stars with x p = 0.1 and —1 ^ a ^ 1. Th e range 
of the symmetry energy is constrained by the requirement that A xv sho uld be invertible (see e.g., Prix et al. 20021 ) as well 
as realistic calculations for neutron star EoS (jLattimer fc Pra kash 2007). For simplicity, we will assume that the symmetry 
energy is constant throughout the star. 



5.2 Model B 



Our second model equation of state is constructed in such a way that we can explore the relevance of the chemical cou- 
pling betweeii_t4ie_two fluid de£rees_of_ft^eedoin_that arises due to composition variation. We consider the simple analytical 
EoS (jPrix fc Rieutordlbooj Undersson et al.ll2002h : 



So = k n pl n + k p pp p 



(49) 



where fc x are constant coefficients. For this EoS, the symmetry energy term vanishes, as A np = 0. However, the polytropic 
indices N x = (7 X — 1) _1 of the neutron and proton fluids can be different, i.e. 7V n 7^ Ap. This enables us to construct stratified 



In compiling Table 1, we noticed an error in the data reported in IPassamonti et al.l (|2009l 1. The data given in the fourth column of Table 
1 of that paper is incorrect. For the sequence of stellar model A, the correct value of the break-up angular velocity is fl^ / \/Gp c = 0.7252 
and the present table gives the correct values for the related quantity Si/Six- We also note a typo in the label of the fifth column in the 
previous paper, where T/\W\ X 10 -2 should be replaced by T/\W\ X 10 2 . 
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Table 2. This table provides the paramete rs for two non-rotating stellar models for the EoS l|49|l . We refer to these models as B NS and 
BO. They correspond to models I and II of iPrix fc Rieutordl (|2002| 1 . respectively. The units of the coefficients fc x are GR\ q p\ ' x , where 
G is the gravitational constant and R eq is the equatorial radius. The proton fraction at the star's centre is x p (0), while the central mass 
density is p c . 



Model 


7n 


7p 


k n 




x p (0) 


Ml (p c R\ q ) 


B NS 


2.0 


2.0 


0.705 


6.347 


0.1 


1.273 


BO 


2.5 


2.1 


0.693 


8.871 


0.1 


1.833 



configurations. To see this, consider the relation between the chemical potential and the mass density 

/ fix 1 " 

px = 



which for a model in /^-equilibrium leads to the following profile for the proton fraction: 

NiV D n-1 



X n 



1 + 



(7nfcn) 

For a vanishing symmetry energy, equations (|36f> — <|3Tf) become 

8P = [(a£- 
S(5 = [(x p - 1) A 

For model B the A 



~N a - 



N„ 



(50) 



(51) 



2x p + 1) A 



1111 H~ 3?p-<4p 



- Zp^pp] Sp + (A pp 



i|,pj pSp + [(x p — 1) Ann + X p App]pSXp , (52) 
+ Ann) S X p ■ (53) 

i xy coefficients are explicitly given by 

Axx = fcx7x(7x-l)Px x " 2 - (54) 
We have constructed two different non-rotating models for this analytic EoS. These two models are labeled B NS and 
B0, and correspond to (after a transformation of units) models I an d II of Prix fc Ricutord (|2002l ). The parameters for 
models B NS and B0 are given in Table [2] Comparing model I and II of IPrix fc Rieutordl (|2002T ) to our numerical models, we 
find an agreement to better than one percent for the dimensionless stellar mass Mj (p c i?g 9 ). 



The non-rotating B 



model represents a non-stratified star, 7 P 

hxx 



7„, with a constant proton fraction given by 



= 0.1. 



k n -f- kp 

This means that, the co- and counter-moving perturbations are decoupled and equations (|52p - (|53[) become: 

kn kr> 



5P = 2- 



-pSp, 



5(3 = 2 (kn + kp) S Xp 



(55) 



(56) 



k n kp 

We use this model to test our evolutions against the frequency domain results of Prix fc Rieutord 1 2002h . 

In addition we consider a rotating sequence, which extends B0 up to the fastest rotating model B12. All our rotating B 
models correspond to the same j x and fe x as the B0 model. Therefore, any rotating B model is stratified with x p (0) 7^ at 
the centre and zero proton fraction at the star's surface. Due to the effect of rotation on the central density and chemical 
potential, the central proton fraction is a; P (0) = 0.1 for the non-rotating model B0 and becomes x p (0) ~ 0.081 for model B12, 
which is near the mass shedding limit. These stratified models are physically interesting, as the relations (|52[1 - (|53[) and the 
perturbation equations couple the co- and counter-moving degrees of freedom and we can study the effect of this coupling on 
the spectrum of stellar oscillations. The main properties of the B models are given in Table [1] 



6 RESULTS 

In this section we present the result from the first ever time-evolution study of perturbed, rapidly rotating, superfluid stars. 
The main aim is to explore how fast rotation, symmetry energy and entrainment affect the non-axisymmetric oscillation 
modes. 



6.1 The evolution code 



The evolution problem for equations (|16[l - (|19p is intrinsically a three dimensional problem. As already mentioned, the Fourier 
decomposition of the azimuthal degree of freedom, identifying specific m-modes, reduces the problem to two spatial dimensions. 
In spherical coordinates, we could then evolve this system of equations on a two-dimensional grid based on the coord inates r and 
8. However, instead of using r we adopt a new radial coordinate x = x(r, 6), fitted to surfaces of constant pressure (jJones et al. 



20021 ; IPassamonti et alj |2009 ) . This leads to an easier implementation of the surface boundary conditions. As we are working 
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Figure 1. This figure shows the frequencies of the axial- led inertial modes and their dependence on the star's spin and the symmetry 
energy a for the sequence of rotating models A with proton fraction x p = 0.1 and entrainment parameter e = 2/3. The mode-frequencies 
and the stellar angular velocity are given in units of y/Gp c , where G is the gravitational constant and p c the central mass density. The 
mode-frequencies are determined in the frame of the rotating star. In the left panel, we show the ordinary r -mode and the superfluid 
r s -modes for different values of a in the range —1/2 $C a sC 4/5. For fast rotating models, Q / \JGp c > 0.25, the r s -mode deviates from 
the simple slow-rotation relation 161> . where 7 S = 1/(1 — e) = 3, because of the symmetry energy. In the right panel we show three 
I = 4, m = 2 superfluid inertial modes for a = —0.5 and a = 0.5. These results also show a clear dependence on the symmetry energy. It is 
worth noticing that, for the non-stratified models A the ordinary inertial modes are equal to the results for single fluid N = 1 polytropes, 
where N is the polytropic index. 



in the time domain, the various mode frequencies are extracted by a Fast Fourier Transformation (FFT) of the time evolved 

perturbation variables. 

(2009) 



The numerical code for superfluid neutron stars extends the single fluid code developed by IPassamonti et al 



In 

fact, with our chosen variables the Euler equation (|16|l and the mass conservation equation (|18|l are identical to the single fluid 
case. We have thus extended the previous code by adding equations (I17|l and (|19|) . together with the appropriate boundary 
conditions. The main elements of the code are the use of a Mac-Cormack algorithm, a second order accurate numerical scheme 
both in space and time, and the implementation of a fourth-order Kreiss-Oliger numerical dissipation, which stabilizes the 
simulations against spur ious high frequency oscil lations. The performance of the final code is practically identical to that of 
the single fluid code, see Passamonti et al. (|2009h for details. 



6.2 Initial Data 



In a time domain study, the initial perturbations can be chosen to excite specific parts of the spectrum. Of course, a strict 
selection of oscillation modes requires the determination of eigenfunctions to be used as initial data. To achieve this, one 
would have to either solve the frequency domain version of the perturbation equation s (I16l) - (|19|) or perform an eigenfunction 
recycling study of the time evolutions ( Stergioulas et alj |2004| ; iDimmelmeier et al.l |2006') . In this work we use neither of 



these strate gies. We are interested in a multi-mode analysis where many oscillation modes are excited in each simulation. As 
discussed by Passamonti et al. ( 20091 ) . this is easily achieved by the use of initial perturbations with an arbitrary radial profile 
and an angular dependence appropriate for the general class of eigenfunction. For the radial part we typically use a Gaussian 
distribution. Meanwhile the angular functions are inspired by slow-rotation results. 

Type I parity perturbations are generally excited with the following mass density and proton fraction perturbations: 



8p — Ao exp 
s Xp = M exp 



r - r 
qR{9) 



( r - r 



Yll (M) , 
Yn (0, cb) , 



(57) 
(58) 



where Ao and A\ are two arbitrary constants that determine the initial values of Sp and 5\ P on the star's surface. The stellar 
radius at polar angle 9 is denoted by R(8), and the parameters ro and q respectively determine the centre of the Gaussian 
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Figur e 2. In this figure we compare our rotating frame frequencies for the r-modes with the slow rotation results of lHaskell et al.l 
(2009), which include O (f! 3 ) corrections. In the left panel, we show the I = m = 2 r s and r -mode frequencies for the rotating models 
A with a = ±1/2 , proton fraction x v = 0.1, and entrainment parameter e = 2/3. Our numerical results are shown as solid lines, while 
the frequencies of lHaskell et alj l|2009h are represented by dashed lines. The right panel provides the relative error between the r s -mode 
frequencies calculated with the two methods, where the dotted lines denote the 2% level. The results show that the calculated frequencies 
agree to better than 2% up to a stellar angular velocity of Q/y/Gp c ~ 0.4. 



profile and its width. The I = m spherical harmonic Yu {0, <j>) approximates the typical angular behaviour of a polar mode in 
a spherical star. For simplicity, all other perturbation variables are set to zero to complete the initial data. 
For Type II parity perturbations we use the following initial data for the vector fields f and D: 

2" 



= p exp 



ro 



D 



(qRW) 



i p exp 



V/ 



/ r - Tp 
\qR(6) 



Yu (0, 



(59) 



(60) 



Here Yjf (6, <j>) is a magnetic spherical harmonic ( Tho me 1980). The remaining perturbations are set to vanish on the initial 
time slice. 



6.3 Inertial modes 



Let us first consider the inertial modes which in a superfluid star split (more or less clearly depending on the equation of 
state) into ordinary and superfluid modes. In the first class, the perturbed fluid elements of the two components oscillate in 
phase, whereas for the superfluid modes they pulsate in counter- phase. As for single fluid st ars, each inertial mode can be 
classified by its parity as an axial- led or polar-led inertial mode ( Lockitch fc Friedman 19991 ). Among the axial-led inertial 
modes, the r-modes form a well defined sub-set. They are the only modes that are purely axial in the slow-rotation limit. In 
stratified neutron stars, only the co-moving r-mode exists. The counter-movin g mode is no longer purely axial, but acquires a 
polar component already at leading order in the rotation l Haskell et al. 2009h . In a non-stratified neutron star, on the other 
hand, the ordinary and superfluid degrees of fre edom are completely decoupled and a purely axial counter-moving r-mode 
exists jAndersson et al.ll2008l ; lHaskell et alJboogh . 



Up to O (fi), the frequencies of the superfluid and ordinary inertial modes are related according to IjPrix et al.1 2004): 

u a ~ 7 e uj , (61) 

where j e = (1 — e) _1 . It makes sense, as a first step, to investigate to what extent the inertial-mode frequencies deviate 
from this simple scaling law in the case of rapid rotation. To do this, we consider model A with fixed entrainment parameter 
e = 2/3 and proton fraction x p = 0.1. For the symmetry energy term we consider four different values a = —1/2, 0, 1/2, 4/5. 
The results are shown in figure [1] In the left panel we show the I = m — 2 ordinary and superfluid r-modes. The ordinary r 
mode is represented by a solid line and the expected frequency (|61(l of the superfluid r s mode is also indicated. Our results 
show that the r s mode frequencies for different values of a agree well with equation (|61[1 roughly up to a stellar rotation 
Q/^Gpc ~ 0.2. For faster rotation, the r s -mode frequency depends strongly on the symmetry energy. 
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Figure 3. This figure illustrates the dependence of the I = m = 2 acoustic modes on the entrainment parameter e for two non-rotating 
stellar models, the non-stratified B NS model (left panel) and the stratified B0 model (right panel). Solid and dashed lines represent 
oscillation modes where neutrons and protons move, respectively, in phase (ordinary modes) and counter- phase (superfluid modes). The 
spectrum of the B0 model shows some avoided crossings. Avoided crossings are expected in a stratified neutron star as the co- and 
counter-moving degrees of freedom are coupled. 



In order to con firm these results, we also considered the superfluid r-modes within the slow-rotation approximation 
20091 ) . The approximate results confirm that, even though one should expect the counter-moving inertial 
modes to approach (|61|) as a — > 1, the relation does not hold perfectly even in the extreme limit. 



( Haskell et al 



Our results show that the r s -mode frequency remains closer to the values expected from equation (|61[) for larger a. Similar 
behaviour is noted for other inertial modes. In the right panel of figure [TJ we show three I — 4, m = 2 axial- led inertial modes 
for a — ±1/2. The dependence on the symmetry energy term is distinguishable beyond fl/^/Gp c ~ 0.25 also for these modes. 
These results can be understood from a local plane- wave analysis, see Appendix [X] for a detailed discussion. 

We can also compare our results to the superfluid r-mode frequencies calculated by Haskell et al. ( 20091 ) . who worked in 
the slow-rotation approximation keeping terms up to O (^ 3 ). The I — m = 2 r s frequency is then given by 

LOs = CO fi + c 2 ft 3 , (62) 

where Co and C2 are constants that depend on the stellar model and the multipole of the modes. The first coefficient has the 
well known analytical expression: 

C0 = 7£ 7(zTT)' (63) 

while C2 is given in closed form by Haskell et al. 1 20091 ). For the I = m = 2 r s -mode of a model A star with e = 2/3 and 
x p = 0.1, we have Co = 2 and the values for C2 given in Table [3] In figure [2] we compare the r s -mode frequencies extracted 



m 

m 



from our evolutions to those determined analytically by Haskell et al.l ((2009J) for different values of the symmetry energy. In 
the left panel, we show how the r a -mode frequency depends on the star's rotation for two selected models with a — ±1/2. The 
right panel shows the relative error between the frequencies obtained with the two methods. The agreement is better than 
three percent up to Q,/y/Gp c « 0.55, which corresponds to a rapidly rotating star with Q/Qk ~ 0.77 (see Table[l}. For faster 
rotation, the errors become larger and reach eleven per cent near the mass shedding limit. In order to improve the accuracy, 
the slow-rotation analysis would need to consider higher order corrections. This may be prohibitively difficult. 

The dependence of the superfluid inertial modes on the proton fraction is very weak for model A neutron stars. That 
this should be expected can be seen from equation (|A9[) . To confirm this we have evolved the fast rotating model A8 with 
x p — 0.01. The frequencies for this model differ by less than one percent from the x p — 0.1 case. 



6.4 Acoustic modes 

Let us now study the effects of entrainment and symmetry energy on the fundamental and pressure modes of a superfluid 
star. The acoustic mode spectrum is characterized by the usual two classes of ordinary (co-moving) and superfluid (counter- 
moving) modes. We will focus on the quadrupole I — m — 2 modes. These are the modes that tend to be the most relevant 
for gravitational- wave studies. 
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Tabl e 3. The c oefficient C2 required in equation Il62jl for the I = m = 2 superfluid r-mode. The values of C2 are given in the units used 
in lHaskell et al, | J2009h . i.e. GM/R 3 . Here we have used the mass M and the radius R for the non-rotating model AO, with e = 2/3 and 
x p = 0.1. 



a C2 



-0.5 -1.052 

0.0 0.452 

0.5 0.956 

0.8 1.123 



First we consider the oscillations of two non-rotating equilibrium configurations, the non-stratified B NS model and the 
stratified B0 model described in Section 15.21 In figure [3] we show the variation of the fundamental quadrupole mode and 
the first two quadrupole pressure modes with the entrainment parameter e. In order to determine the co- and counter- 
moving character of a mode, we have first reconstructed the time variation of the component velocities v x from the primary 
dynamical variables f and D. Post processing the numerical evolution dat a, we have determined the veloci ty components using 
the eigenfunction extraction code developed by Stergioulas et al. (2004) and Dimmelmeier et al. ( 20061 ) . For any mode, we 



have then compared the velocity eigenfunctions of the two fluid components and determined whether they oscillate in phase 
or counter-phase. In figure [3l the solid lines represent modes that oscillate in phase, while dashed lines correspond to modes 
that pulsate in counter-phase. In absence of composition gradients (model B NS in the left panel of figure [3}, the co-moving 
and counter-moving degrees of freedom are completely decoupled and the spectrum does not exhibit any interaction between 
ordinary and superfluid modes. They are actually related by the simple expression, cf. (|A9[) . 



lj s ~ /ifiiJo • (64) 

The spectrum of stratified superfluid stars is more interesting. From the right panel of figure [3] we note that the coupling 
between "ordinary" and "superfluid" perturbations generates avoided crossings, where the oscillation phase of the mode 
changes at (roughly) e = 0. This behaviour is more evident in the ordinary and superfluid pressure modes. From our data 
there does not appear to be an avoided crossing for the fundamental mode. However, the region where this crossing should 
take place is difficult to resolve with our evolutions. We have tried different initial data sets but only managed to find a singl e 
f-mode at e = 0. Nonradial oscillations of the B NS and B0 models have already been studied by Prix fc Rieutordl ( 20021 ) . 



although not within the Cowling approximation. A direct comparison with their results is not possible since the Cowling 
approximation can introduce a 15-20% error in the f-mode frequencies. However, the qualitative behaviour of the acoustic 
spectrum in the two studies is clearly similar. 

In order to investigate the effect of entrainment on the rotational splitting of the non-axisymmetric acoustic modes we 
consider a sequence of A models with a — and models B. In figure [4] we show results corresponding to the I — m — 2 acoustic 
modes in the e = 0.5 case. For the generic initial data given by equations (|57p - (|58[) many oscillation modes are excited in 
the numerical simulations. By studying their eigenfunctions we can track individual modes from the non-rotating model up 
to the mass shedding limit. An example is illustrated in the left panel of figure [4] where the rotating frame frequencies of 
the superfluid f-mode and the first two p-modes are shown for the A models. The counter-moving f s mode exhibits a larger 
rotational splitting compared to the two pressure modes. In the right panel of figure [4] we compare instead the normalized 
frequencies of the superfluid f-mode, namely uj/uj^ 11 , for models A and B. The horizontal-axis of this figure represents the 
stellar rotation divided by the maximal angular velocity Qk ■ The results in the figure show that, even though the two classes 
of models have the same value of the entrainment parameter, the non-axisymmetric splitting of the f s mode is quite different in 
the two cases. This is particularly clear in the rapid rotation regime. Meanwhile, the two f s modes have very similar frequencies 
for Q < 0.25 Q.K- 
it is instructive to explore the slow-rotation regime to see if there is a simple relation between the f s mode frequencies 
and the entrainment. To this end, we determine the f s mode frequencies for models A and B with —0.7 ^ e ^ 0.7. In the left 
panel of figure the variation of the f s mode frequency with e is shown for models A. It is useful to recall that these models 
have a = by construction. All modes have been normalized to the f-mode frequency of the non-rotating star. The rotational 
splitting of the modes strongly depends on the parameter e. In order to quantify the effect, we approximate the dimensionless 
mode- frequency lu b = ui B /^/Gp c (see Appendix |A")) by a second order polynomial in Cl — Q/^/Gp c , 

uis = ^ K + ci (e, a, m) Cl + c 2 (e, a, m)Cl 2 + O (&J > (65) 

where Lj^ n = uj^ r j\JGp c is the superfluid mode frequency of the non-rotating star, while ci and C2 are two parameters that 
we fit to the numerical data. In general, these parameters are functions of the entrainment, the symmetry energy and the 
multipole m of the mode. Since this is inherently a slow-rotation approximation, we will only consider rotating models with 
O <^ 0.375 SI a'- An accurate description of the spectrum for more rapidly rotating models would require higher order fitting 
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Figure 4. This figure shows the rotational splitting of the I = m = 2 non-axisymmctric modes (in the reference frame of the rotating 
star) for the two sequences of rotating models A and B. In the left panel, we show the dimensionless frequencies of the superfluid f-modc 
and the first two counter-moving p-modes for models A with cr = and e = 0.5. In the right panel, we compare the superfluid f-modcs 
of these A models with the models B having E = 0.5. In this panel, the star's angular velocity is normalized (the horizontal axis) with 
the mass-shedding limit Q/fiji, while the mode-frequency is divided by the superfluid f-mode frequency of the non-rotating star oj^ 13 
(the vertical axis). 



functions. In the right panel of figure \5\ we show the dependence of uj^ r and ci on 

7e = (l-e)- 1 , (66) 

for Models A. We have denned c[ and for retrograde and prograde modes, respectively. If we assume that a perturbation 
variable is proportional to e H l,ut + m W f modes with positive (negative) m move retrograde (prograde) with respect to the stellar 
rotation. 

For non-stratified and non-rotating stellar models, equation (IA9I) suggests that the superfluid and ordinary f-mode 
frequencies are related according to 

. NR , i \ - NS /a n\ 

^ s = V7e Xto Xp) Wo , (67) 
where \ is a function of a and x p that, in general, depends on the EoS. For models A, this function is given explicitly by 

X» S £±gj^ . (68) 

1 - (1 + <T)Xp 



Since we consider a sequence of models with a — it follows that x — !• We test the accuracy of equation (|67[) against the 
f-mode frequencies for models A shown in figure [5] Fitting our numerical data to a function of form 

ajl , (69) 



we obtain the values for a and 6 listed in Table [4] These results are in very good agreement with the analytical formula (|67l) . 
In fact, the ordinary f-mode frequency determined from our code is 0> o = 1.867, which agrees to better than 0.06% with the 
fitted value. Equation (|67|l also provides an accurate representation for the f-mode of the B models. In this case, the numerical 
evolutions lead to Q — 2.1598, which is within the error bar of the fitted result, see Table 2] It should be noted that we do 
not have a simple analytic expression for \ m t ne case of the B models. Our results suggest that x is only weakly dependent 
on the entrainment also for this equation of state. 

Let us now consider the rotational corrections to the counter-moving f-mode frequency, cf. equation (|65|l . We find that 
we can still use the fitting function (|69p for the parameter ci. This leads to the results given in Table [4] These results suggest 
that the rotational splitting parameter ci depends linearly on the entrainment parameter 7 e for both models A and B. This 
is also clear from the results in Figure [5] To describe the quantities c£ and cf, we instead used the following fitting function: 

a + b ^ . (70) 

£NR W 

The results for a and b are given in Table [4] Despite the fact that models B are stratified, the f s mode seems to maintain the 
same dependence on the entrainment as for the non-stratified models A in the slow rotation regime. 

Building on this, let us try to understand the dependence on the symmetry energy. To do this we consider the sequence 
of rotating stellar models A and take the symmetry energy in the range — 1/2 ^ cr ^ 1/2. Meanwhile the proton fraction and 
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Figure 5. The effect of the entrainment parameter e on the rotational splitting is shown in this figure. For the sequence of stellar models 
A with cr = 0, the left panel displays the rotating frame frequencies of the I = m = 2 superfluid modes for several values of e. The 
mode- frequencies are normalized with the superfluid f-mode frequencies of the non-rotating star In the right panel, we show how 

the dimensionless frequency uj^ 1 ^' = uj B /y/Gp c and the two splitting parameters c\ and Cj of equation H65II depend on the entrainment 
parameter j s = (1 — e) — . The superscripts r and p for the parameters cj and Cj denote retrograde and prograde modes, respectively. 
The quantities Cj and c? 



the entrainment are kept constant at x p = 0.1 and e = 2/3, respectively. The f s mode dependence on the stellar spin (for 
some of these models) is shown in figure [6] In the figure, we have normalized the mode-frequencies to the f-mode frequency 
of the non-rotating model. The results show that the rotational splitting of the I = m = 2 superfluid f-mode decreases for 
larger values of the symmetry energy. In the slow rotation regime, we can study this effect by fitting the numerical data using 
equation (|65|) . Adding the dependence on the symmetry energy to our previous results, we now use the fitting function 

= x(°"i x p) &o S + 7e ci (a, m) Cl + c 2 (s, a, m) Cl 2 + O [Cl 3 j , (71) 

where £2 is in general different from C2, and we have assumed that c\ is a function of only a and m. First we consider the 
non-rotating case, using the ordinary frequency of the f-mode d)^ s as a fitting parameter. The fitted result agrees with the 
numerical value determined from the simulations to better than 0.04%. Empirically we find that the parameters £1 can be 
approximated by the function: 

(72 ) 

The determined values for a and b for model A are given in Table [5] We find that the parameter S2 assumes less regular 
values. They can be fitted with a quadratic polynomial in a, but the fits are not very accurate. Hence, we do not provide any 
of those results here. However, as our main aim was to provide an approximate description of the rotational splitting of the 
quadrupole superfluid f-mode, reliable values for the coefficient £1 should be sufficient. 



7 CONCLUSIONS 

In this paper we have considered, for the first time, the oscillations of a superfluid neutron star as an initial-value problem. Using 
time evolutions of the relevant linearised equations we studied non-axisymmetric oscillations of rapidly rotating superfluid 
neutron stars. We considered perturbations of axisymmetric background configurations in Newtonian gravity and accounted 
for the presence of superfluid components via the standard two-fluid model. Within the Cowling approximation, we were able 
to carry out evolutions for uniformly rotating stars up to the mass-shedding limit. Our results represent the first detailed 
analysis of superfluid neutron star oscillations in the fast rotation regime, where the star is significantly deformed by the 
centrifugal force. 

For simplicity, we focused on background models such that the two fluids (superfluid neutrons and protons) co-rotate, are 
in /3-equilibrium and coexist in all the volume of the star. Two different analytical model equations of state were considered. 
The models were chosen to represent relatively simple generalizations of single fluid, polytropic stars. We investigated the 
effects of entrainment, rotation and symmetry energy on various non-radial oscillation modes of these models. Our results 
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Table 4. This table provides the values of the fitting parameters a ± Aa and b ± Af> from equations J69H and (1 70 6 for the superfluid 
f-mode of the non-rotating star, i.e. and the splitting coefficients c r k and (? h with k = 1, 2 for the retro- and pro-grade f a modes. In 

the fits we have used the stellar models A and B with Q < 0.375 Q.k ■ 





Models 


a 


Aa 




b 


A6 




^s NR 


A 


1.8686 


2.1 x 10" 


-4 


0.5000 


1.5 x 10" 


4 




A 
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-3 




B 


-0.9837 


3.8 x 10- 
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1.0059 
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-3 


r p 
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-0.1461 


1.2 x 10- 


-2 


0.1106 


1.1 x 10- 
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show that entrainment and symmetry energy can have a significant effect on the rotational splitting of non-axisymmetric 
modes. In particular, the symmetry energy modifies the inertial mode frequencies considerably in the regime of fast rotation. 

The perturbative time-evolution framework provides a useful tool that should allow us to consider more realistic (and 
by necessity complicated) neutron star models in the future. The clear advantage over frequency-domain studies is that it 
is straightforward to study oscillations corresponding to eigenfunctions with a complex set of rotational couplings. This is 
particularly useful in the rapid rotation regime. The obvious downside is that time evolutions can never provide the "complete" 



Table 5. This table provides the values of the fitting parameters a ± Aa and b ± Af> from equation H72I I. They describe the splitting 
functions and Cj for the I = m = 2 superfluid f-mode. The labels r and p denote retro- and pro-grade f s modes, respectively. In the 
fits we have used results for models A with Q < 0.375 Ok- 
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1.1 x 10- 3 
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1.5 x 10~ 2 


5? A 


-0.3045 


1.5 x 10" 3 


-1.024 


2.0 x 10" 2 
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mode spectrum of the star. Initial data has to be chosen in such a way that the oscillations of interest are excited at a significant 
level. It is difficult to, without prior knowledge, find initial data that excites only a few modes. In many cases this is, however, 
less relevant. The main question is if one can extract accurate information regarding the nature of the star's oscillations from 
the numerical data. The results we have presented demonstrate that this is, undoubtedly, the case. Hence, it is relevant to 
develop the perturbative evolution framework further. We are currently considering more general background models, with a 
relative velocity between the two fluid components. We are also adding the dissipative coupling associated with mutual friction 
to the code. Once these features are incorporated we will be able to consider (obvio usly still at a basic level) the dynamics 
associated with the superfluid two-stream i nstability ( Andersson et al. 20031 . 2004 ) and the possible relation with pulsar 
glitches (see Glampedakis fc Anderssonll2008l . for a recent discussion). This is a very exciting prospect. Looking further ahead, 
we would like to add layering to the stellar model by introducing both an elastic crust and distinct superfluid/normal regions. 
There are challenges associated with these aspects, but there is no reason why these developments should be prohibitively 
difficult. 



ACKNOWLEDGEMENTS 

This work was supported by STFC through grant number PP/E00I025/I. 



APPENDIX A: LOCAL ANALYSIS 



In this Appendix we carry out a local analysis of the perturbation equations (|16p ~ (|19|) . The aim is to obtain understand the 
superfluid oscillation spectrum and its dependence on parameters like the proton fraction, the entrainment and the symmetry 
energy. Since we are only interested in a qualitative picture we focus on the non-stratified N = 1 polytropic model A (see 
Section r5.1|) . For this sequence of models, the ordinary and superfluid perturbation variables are decoupled. We c onsider only 



the " superfluid" perturbations as the dispersion relations for the ordinary modes are well known, see for example (jUnno et al 
19891 ). For model A the superfluid perturbation equations (|17[) and (|19[) can be written 



77 1 d t D = 
dtSxp = 
where we have defined 



a 2 V<5x P — 2J7 x D . 
V ■ D - f ■ Vx B , 



(1 



c 2 s (1 + a) (1 - x p ) 

[l-(l + ff)Xp] 



(Al) 
(A2) 

(A3) 
(A4) 



and where the speed of sound for an N = 1 polytrope is given by c 2 s = 2Kp. Now we assume that the perturbation variables 
behave as plane waves, i.e. we introduce an e l '" t_k r ' dependence for all perturbations into equations (|A1[) - ()A2[1 . Here u> and 
k are the frequency and wave vector, respectively. The characteristic polynomial of the resulting equations is then given by 

, 2 



/ -2 . . 2A2\ ~2 . ~2 3 

( V 7e + 47 e 12 I u + rj 7 e 



(2fl ■ kj' 



= 0, 



where k is the unit wave vector and we have defined the following dimensionless variables: 
Q 



n = 



ak 

7Wc 



(A5) 

(A6) 
(A7) 
(A8) 



The quantities M and R denote the mass and radius of a non-rotating stellar model, respectively. In the slow-rotation 
approximation, we can assume that f2 << i). Equation (|A5|I then have the following solutions [up to O (fi 3 )] 

3/2 

'fi 2 -fn-kT , (A9) 



L02 



V7eV + 

7 £ 2fi ■ k 



r, 3 

2 7 e 



V 

4 7 2 



Ok 



n • k ~ 7 e 2fi • k \ 1 



2 7e 



ftk 



(A10) 



As an estimate we consider a non-rotating background model with an = 1 polytropic EoS. Using the definitions (|A4[I 
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T 



superfluid inertial modes 




Figure Al. This figure shows the inertial mode frequencies estimated from the dispersion relation | |A9| I, We consider the four values of 
the symmetry energy term given in the legend. In equation |)A9|I . the correction term at C(f2 3 ) leads to a dependence on a which is very 
similar to the behaviour seen for the numerically determined mode-frequencies, see figure [T] in the main text. 



and (IA8I) we have 



48 (1 + ^(1-^) (RV 
77 " n [l-(l + a)x p ] [x] ' (A11) 
where A is the wavelength (related to the the wave vector according to k = 2n/X). In addition, we have replaced the speed of 
sound with its average value for an N = 1 polytrope: 

<c2> = 4^. (A12) 
■k z R 

If we now assume that the wavelength of the mode is of the same order as the stellar radius, A = R, equation (jAllfl for 
i) becomes 

2 ^48(l + a)(l-x P ) 

' 7T [1 - (1 + Cf) X P ] K ' 

To get a qualitative picture, we consider parameter values e = 2/3 and x p = 0.1 and specify the angle between f2 and k to 
be 8 = 7r/4. In figure |A"T1 we show the positive frequency a; s of equation (|A9f) for different values of a. The O (^ 3 ) correction 
term leads to a behaviour that resembles the global mode results shown in figure [1] in the main text. 
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